A few main packages are dedicated to spatial data handling like import/export, display, geoprocessing, spatial ststistics.
rgdal: interface between R and GDAL (Geospatial Data Abstraction Library) and PROJ4 libraries: raster / vector geospatial data formats and coordinate transformation.
sp: classes and methods for spatial data in R.
rgeos: interface between R and GEOS (Geometry Engine - Open Source) library: area, perimeter, distances, dissolve, buffer, overlap, union, contains…
These packages are still widely used.
sf Website: Simple Features for R
First release: October 20, 2016
sp, rgeos and rgdal functionnalities in one package.
Easier data handling, simpler objects.
Tidy data: compatibility with the pipe synthax and tidyverse operators.
Main author and maintainer: Edzer Pebesma (also sp author)
sf objects data structure:
sfReading layer `martinique' from data source `C:\Users\kim.antunez\Desktop\Florence_R\flo\data\mtq\martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID): 32620
proj4string: +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add = TRUE, cex = 1.2, col = "red", pch = 20)Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.000 35297.56 3091.501 12131.617 17136.310
[2,] 35297.557 0.00 38332.602 25518.913 18605.249
[3,] 3091.501 38332.60 0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702 0.000 7177.011
[5,] 17136.310 18605.25 20226.198 7177.011 0.000
Simple union:
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2, border = "red")Aggregation according to a grouping variable:
library(dplyr)
mtq_u2 <- mtq %>% group_by(STATUT) %>% summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u2), add = T, lwd = 2, border = "red", col = NA)mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2)
plot(st_geometry(mtq_b), add = T, lwd = 2, border = "red")m <- rbind(c(700015, 1624212), c(700015, 1641586), c(719127, 1641586), c(719127,
1624212), c(700015, 1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border = "red", lwd = 2, add = T)mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = T)google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)
mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join = st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col = "lightblue")Several solutions are available:
ggplot2 users can have a look to ggplot2 mapping features (geom_sf) that can mix nicely with ggspatial.tmapcartography is based on base graphics and allow most of basic and advanced cartographic representations. Full disclosure: the speaker is the maintainer of cartography.Here we will focus on cartography and do small examples with ggplot2.
cartographylibrary(sf)
# Import geo layer
sm <- st_read(dsn = "../data/rhone/seine_maritime.geojson", stringsAsFactors = F,
quiet = TRUE)
dep <- st_read(dsn = "../data/rhone/dep.geojson", stringsAsFactors = F, quiet = TRUE)
# change projection (lambert93)
sm <- st_transform(sm, 2154)
dep <- st_transform(dep, 2154)
# Import dataset
csp <- read.csv("../data/rhone/data_seine_maritime.csv")
# merge geolayer and dataset
sm <- merge(sm, csp, by = "INSEE_COM", all.x = TRUE)library(cartography)
plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")# Custom map of active population
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
bb <- st_bbox(sm)
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
plot(st_geometry(sm), col = "cornsilk2", border = NA, lwd = 0.5, add = T)
propSymbolsLayer(sm, var = "act", col = "darkblue", inches = 0.5, border = "white",
lwd = 0.7, symbols = "square", legend.style = "e", legend.pos = "topleft",
legend.title.txt = "Labor Force\n(2014)", legend.values.rnd = 0)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
typoLayer(sm, var = "cat", border = "ivory", lwd = 0.5, legend.values.order = c("agr",
"art", "cad", "int", "emp", "ouv"), col = c("#e3b4a2", "#a2d5d6", "#debbd4",
"#b5dab6", "#afc2e3", "#e9e2c1"), add = TRUE, legend.pos = "n")
legendTypo(title.txt = "Dominant socio-professional Category", col = c("#e3b4a2",
"#a2d5d6", "#debbd4", "#b5dab6", "#afc2e3", "#e9e2c1"), categ = c("Agriculteurs",
"Artisans", "Cadres", "Prof. Inter.", "Employés", "Ouvriers"), nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm, var = "pcad", breaks = bks, col = cols, border = "grey80", lwd = 0.4,
legend.pos = "topleft", legend.title.txt = "Share of managers (%)", add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers (%)",
add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)These maps may allow to hide or, at least, diminish the MAUP.
3904.689 m
grid <- getGridLayer(x = sm, cellsize = 3900 * 3900, type = "regular", var = c("cad",
"act"))
grid$pcad <- 100 * grid$cad/grid$act
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(grid, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers\n(en %)",
add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)grid$cad100 <- grid$cad * 100
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", typefct = "exponential",
span = 4000, beta = 2, breaks = bks, col = cols, legend.pos = "topleft",
mask = st_buffer(sm, 0), legend.title.txt = "Share of managers* (%)", border = "grey90",
lwd = 0.2, add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)
text(x = 518000, y = 6921123, font = 3, cex = 0.8, labels = "*Potential smoothing\n exponential function\n span = 2.5km, beta = 2")ggplot2library(ggplot2)
ggplot() +
geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
geom_sf(data = sm %>% st_centroid(),
aes(size= cad), colour="#E84923CC", show.legend = 'point') +
scale_size(name = "Active population",
breaks = c(100,1000,5000),
range = c(0,20)) +
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(sm)[c(1,3)],
ylim = st_bbox(sm)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(title = "Active population",
caption = "Insee, 2018\nKim & Tim, 2018")ggplot() +
geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
geom_sf(data = sm, aes(fill = pcad), colour = "grey80") +
scale_fill_gradientn(name = "Share of managers (%)",
colours = carto.pal("green.pal", 3,"wine.pal",3),
values=bks/max(bks))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(sm)[c(1,3)],
ylim = st_bbox(sm)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(title = "Managers",
caption = "Insee, 2018\nKim & Tim, 2018")CRAN Task View: Analysis of Spatial Data
See file data_prep.R for data extraction.
# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)# Load data
adm <- readRDS("../data/iris_p13.rds")
feat_sir <- readRDS("../data/sir_p13.rds")
feat_osm <- readRDS("../data/osm_p13.rds")
# Define margins
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
# Plot the maps
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "SIR", sources = paste0(nrow(feat_sir), " restaurants"),
author = "Kim & Tim, 2018", scale = NULL, frame = FALSE, tabtitle = TRUE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, bty = "n",
cex = 0.7, col = "red")
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "OSM", sources = paste0(nrow(feat_osm), " restaurants"),
author = "Kim & Tim, 2018", scale = 0.5, north = T, frame = FALSE, tabtitle = TRUE)Explore datasets interactively
The default is not bad:
mapview allows highly customizable maps:
mapview(feat_sir, map.types = "OpenStreetMap", col.regions = "#940000", label = paste(feat_sir$L2_NORMALISEE,
feat_sir$NOMEN_LONG, sep = " - "), color = "white", legend = TRUE, layer.name = "SIRENE",
homebutton = FALSE, lwd = 0.5) + mapview(feat_osm, col.regions = "#000094",
color = "white", legend = TRUE, label = feat_osm$name, lwd = 0.5, layer.name = "OSM",
homebutton = FALSE)These maps, as appealing as they seem to be, are not really suitable for displaying geostatistical information.
Nonetheless they can be really useful for exploratory data analysis.
# Intersect adm and feat
inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)
# Count points in polygons
adm <- st_sf(adm[, 1:2, drop = T], n_osm = sapply(X = inter_osm, FUN = length),
n_sir = sapply(X = inter_sir, FUN = length), geometry = st_geometry(adm))
# Display the map
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, tabtitle = TRUE,
frame = FALSE)
# Display the map
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T,
tabtitle = TRUE, frame = FALSE)Create a small function that builds a regular grid and counts points in each cell.
pt_in_grid <- function(feat, adm, cellsize = 50) {
# Create a regular grid (adm bbox)
grid <- st_make_grid(x = adm, cellsize = cellsize, what = "polygons")
# Keep only cells that intersect adm
. <- st_intersects(grid, adm)
grid <- grid[sapply(X = ., FUN = length) > 0]
# Count pts in grid
. <- st_intersects(grid, feat)
grid <- st_sf(n = sapply(X = ., FUN = length), grid)
return(grid)
}Create a function that plots a map
plot_grid <- function(grid, adm, title, bks, col) {
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5)
choroLayer(grid, var = "n", border = NA, add = TRUE, breaks = bks, col = cols,
legend.pos = "n")
layoutLayer(title = title, scale = NULL, frame = FALSE, tabtitle = TRUE,
author = "", sources = "")
}Build grids and display maps
# plot 6 maps on a single figure
par(mar = c(0, 0, 1.8, 0), mfrow = c(2, 3), bg = "#FBEDDA", ps = 16)
# defines a unique set of breaks for all maps
bks <- seq(1, 22, 3)
cols <- carto.pal("sand.pal", length(bks) - 1)
# SIR maps 50 m cells
grid <- pt_in_grid(feat_sir, adm, 50)
plot_grid(grid = grid, adm = adm, title = "SIR - 50 m cells", col = cols, bks = bks)
# 100 m cells
grid <- pt_in_grid(feat_sir, adm, 100)
plot_grid(grid = grid, adm = adm, title = "SIR - 100 m cells", col = cols, bks = bks)
# 150 m cells
grid <- pt_in_grid(feat_sir, adm, 150)
plot_grid(grid = grid, adm = adm, title = "SIR - 150 m cells", col = cols, bks = bks)
# Add a north arrow
north()
# OSM maps 50 m cells
grid <- pt_in_grid(feat_osm, adm, 50)
plot_grid(grid = grid, adm = adm, title = "OSM - 50 m cells", col = cols, bks = bks)
# Add sources
mtext(text = "Map data © OpenStreetMap contributors, under CC BY SA.\nSIRENE - 2018, INSEE - 2018",
side = 1, line = -1, cex = 0.5, adj = 0)
# 100 m cells
grid <- pt_in_grid(feat_osm, adm, 100)
plot_grid(grid = grid, adm = adm, title = "OSM - 100 m cells", col = cols, bks = bks)
# 150 cells
grid <- pt_in_grid(feat_osm, adm, 150)
plot_grid(grid = grid, adm = adm, title = "OSM - 150 m cells", col = cols, bks = bks)
# Add legend
legendChoro(pos = "topright", title.cex = 0.7, values.cex = 0.5, title.txt = "N. restaurants",
breaks = bks, nodata = FALSE, values.rnd = 0, col = cols)
# Add scalebar
barscale(1)Define a function that computes smoothed restaurant density (KDE method).
compute_kde <- function(feat, adm, title, sigma = 100, res = 50) {
# Define an observation window
w <- as.owin(as(adm, "Spatial"))
# sf to coords
pts <- st_coordinates(feat)
# Coords to ppp
p <- ppp(pts[, 1], pts[, 2], window = w)
# Compute KDE
dens <- density.ppp(p, sigma = sigma, eps = res)
# Image to raster (+ proj & km2)
result <- raster(dens, crs = st_crs(adm)[[2]]) * 1e+06
return(result)
}Define a function that maps the results
plot_kde <- function(x, adm, method = "equal", nclass = 12, title = "", author = "",
sources = "", scale = NULL, legend.title = "Density\n(n/km2)", legend.pos = "topright",
values.rnd = 0) {
# compute breaks
bks <- unique(getBreaks(values(x), nclass = nclass, method = method))
# Color ramp
cols <- (mapview::mapviewGetOption("raster.palette"))(length(bks) - 1)
# Plot the map
plot(st_geometry(adm), col = NA, border = NA, main = "", bg = "#FBEDDA")
plot(x, breaks = bks, col = cols, add = T, legend = F)
legendChoro(pos = legend.pos, title.txt = legend.title, values.cex = 0.5,
breaks = bks, nodata = FALSE, values.rnd = values.rnd, col = cols)
layoutLayer(title = title, scale = scale, tabtitle = TRUE, frame = FALSE,
author = author, sources = sources)
}# Plot 6 maps on a single figure
par(mar = c(1, 0, 1.8, 0), mfrow = c(2, 3), bg = "#FBEDDA", ps = 16)
x <- compute_kde(feat = feat_sir, adm = adm, sigma = 50, res = 20)
plot_kde(x, adm, title = "SIR - Sigma = 50", legend.title = "Restaurants\n/km2")
x <- compute_kde(feat = feat_sir, adm = adm, sigma = 100, res = 20)
plot_kde(x, adm, title = "SIR - Sigma = 100", legend.title = "Restaurants\n/km2")
x <- compute_kde(feat = feat_sir, adm = adm, sigma = 250, res = 20)
plot_kde(x, adm, title = "SIR - Sigma = 250", legend.title = "Restaurants\n/km2")
x <- compute_kde(feat = feat_osm, adm = adm, sigma = 50, res = 20)
plot_kde(x, adm, title = "OSM - Sigma = 50", legend.title = "Restaurants\n/km2")
x <- compute_kde(feat = feat_osm, adm = adm, sigma = 100, res = 20)
plot_kde(x, adm, title = "OSM - Sigma = 100", legend.title = "Restaurants\n/km2")
x <- compute_kde(feat = feat_osm, adm = adm, sigma = 250, res = 20)
plot_kde(x, adm, title = "SIR - Sigma = 250", legend.title = "Restaurants\n/km2")reproducibility
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] mapview_2.3.0 maptools_0.9-2 spatstat_1.56-1
[4] rpart_4.1-10 nlme_3.1-128 spatstat.data_1.3-1
[7] ggplot2_3.0.0 SpatialPosition_1.2.0 raster_2.6-7
[10] sp_1.2-5 cartogram_0.1.0 cartography_2.1.1
[13] dplyr_0.7.4 sf_0.6-3 rmdformats_0.3.3
[16] knitr_1.20
loaded via a namespace (and not attached):
[1] viridisLite_0.3.0 foreach_1.4.3 R.utils_2.5.0
[4] shiny_1.0.5 assertthat_0.2.0 highr_0.6
[7] stats4_3.3.2 yaml_2.1.18 pillar_1.1.0
[10] lattice_0.20-34 glue_1.2.0 digest_0.6.15
[13] polyclip_1.9-1 colorspace_1.3-2 htmltools_0.3.6
[16] httpuv_1.3.5 Matrix_1.2-7.1 R.oo_1.21.0
[19] plyr_1.8.4 pkgconfig_2.0.1 questionr_0.6.2
[22] bookdown_0.5 xtable_1.8-2 webshot_0.5.0
[25] scales_0.5.0.9000 tensor_1.5 satellite_1.0.1
[28] spatstat.utils_1.9-0 tibble_1.4.2 mgcv_1.8-15
[31] withr_2.1.2 lazyeval_0.2.1 gdalUtils_2.0.1.7
[34] magrittr_1.5 mime_0.5 deldir_0.1-14
[37] evaluate_0.10.1 R.methodsS3_1.7.1 foreign_0.8-67
[40] class_7.3-14 tools_3.3.2 formatR_1.5
[43] stringr_1.3.1 munsell_0.4.3 bindrcpp_0.2
[46] e1071_1.6-8 rlang_0.2.1 classInt_0.1-24
[49] units_0.6-0 grid_3.3.2 iterators_1.0.8
[52] rstudioapi_0.7 htmlwidgets_1.2 goftest_1.1-1
[55] crosstalk_1.0.0 miniUI_0.1.1 labeling_0.3
[58] base64enc_0.1-3 rmarkdown_1.10.12 gtable_0.2.0
[61] codetools_0.2-15 abind_1.4-5 DBI_1.0.0
[64] R6_2.2.2 rgdal_1.2-11 rgeos_0.3-24
[67] bindr_0.1 stringi_1.1.7 htmldeps_0.1.1
[70] Rcpp_0.12.16 png_0.1-7 leaflet_1.1.0